Proyecto - Pensiones I

Librerías e importaciones

source("cod/setup.R")
source("cod/load_database.R")

Gráficos de CDR

source("cod/graficos.R")

Preparación Balance

Población general

cotizantes <- BD_Cotizantes %>% select(Fec.Nac, Sexo)
names(cotizantes) <- c('fecha_nac', 'sexo')
pensionados <- BD_Pensionados %>% select(FEC_NAC, SEXO)
names(pensionados) <- c('fecha_nac', 'sexo')
pensionados <- pensionados  %>%
  mutate(sexo = ifelse(sexo == "M", 1, 2))
poblacion <- rbind(cotizantes, pensionados) 
poblacion <- poblacion %>% mutate(
  edad = 2024 - year(poblacion$fecha_nac)
)
cuentas <- poblacion %>% count(sexo,edad)
unicos <- cuentas %>% select(sexo, edad)

pensionados <- pensionados %>% mutate(
  edad = 2024 - year(pensionados$fecha_nac),
  meses = 12 - month(pensionados$fecha_nac)
)
cotizantes <- cotizantes %>% mutate(
  edad = 2024 - year(cotizantes$fecha_nac),
  meses = 12 - month(cotizantes$fecha_nac)
)
cotizaciones <- BD_Cotizantes %>% select(-1,-2,-3)
pensionados <- cbind(pensionados, BD_Pensionados) 
pensionados <- pensionados %>% select(-c(4,8))
pensionados <- pensionados %>% mutate(
  edad = 2024 - year(pensionados$fecha_nac)
)

niveles <- 110.39017/IPC$Nivel[229:588]

Funciones alternas

source("cod/curva_dens.R")
source("cod/curva_salarial.R")
source("cod/funciones_alternas.R")

Curva salarial

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Propiedades de cotizaciones

cotizantes <- cbind(cotizantes, cot_prop())

Curva de densidades de cotizacion

curv_dens <- curva_densidad$media
df_dens <- data.frame(
  Periodo = curva_densidad$edad,
  Densidad = curv_dens
)

fig <- ggplot(df_dens, aes(x = Periodo, y = Densidad)) +
  geom_line(color = "#009E73", size = 1.2) +
  geom_point(color = "#009E73", size = 2) +
  labs(
    title = "Curva de densidad de cotizaciones",
    x = "Período",
    y = "Densidad"
  ) +
  theme_minimal() 
fig %>% ggplotly()

# ggsave("res/curva_dens.pdf", fig, width = 8, height = 5)

Obtención de matrices de probabilidad de muerte

find_qx <- function(x, sexo, año){
  return(SUPEN[which(SUPEN$sex == sexo & SUPEN$year == año & SUPEN$edad == x),6])
}
año <- 2025
system.time(
all_qx <- list(
  sapply(0:95, function(y) sapply(0:115, function(x) find_qx(x, 1, año+y))),
  sapply(0:95, function(y) sapply(0:115, function(x) find_qx(x, 2, año+y))))
)

##    user  system elapsed 
##   10.81    1.76   13.59

Precálculo de las pensiones

sal_pen <- function(){
  pensiones <- list()
  for(id in 1:5196){
    x <- cotizantes$edad[id]
    cuota <- cotizantes$cuotas[id]
    sal_prom <- cotizantes$sal_prom[id]
    cuotas_past <- cotizaciones[id,]
    n_cot <- c(0,cumsum(curv_dens[(x-19):96]))+cuota
    salarios <- cumprod(curv_sal[(x-19):96])*sal_prom
    enteros <- round(n_cot - n_cot[1])
    cantidades <- enteros[-1] - enteros[-length(enteros)]
    sal_pen <- rep(0, 115-x)
    sal_pen[1] <- salarios[1]
    cuotas <- unlist(cuotas_past*niveles)
    cuotas <- cuotas[cuotas>10000]
    if(sal_prom == 0){ # supone que ya no va a trabajar?
      pensiones[[id]] <- rep(mean(head(sort(
        cuotas[sum(1 - cuotas < 5e6)/enteros[i+1] > 0.5 | cuotas < 5e6],
        decreasing = T), 300)), 115-x)
      next
    }
    for(i in 1:(114-x)){
      cuotas <- c(cuotas, rep(salarios[i+1], cantidades[i]))
      calc_pen <- cuotas[sum(1 - cuotas < 5e6)/enteros[i+1] > 0.5 | cuotas < 5e6]
      sal_pen[i+1] <- mean(head(sort(calc_pen, decreasing = T), 300))
    }
    porc <- c(0.525,0.51, 0.494, 0.478, 0.462, 0.446, 0.43)
    montos <- c(2,3,4,5,6,8) * 367108.55
    indices <- findInterval(sal_pen, montos) + 1
    sal_pen <- sal_pen * porc[indices]
    pensiones[[id]] <- sal_pen
  }

  return(pensiones)
}

system.time(pensiones <- sal_pen())

##    user  system elapsed 
##   71.97    1.03   74.03

Obtención del interés

inf_max <- IPC$Nivel[564]/IPC$Nivel[552] - 1
(int_real <- mean(int) - inf_max)

## [1] 0.07793612

Primer modelo determinístico

Función del modelo

source("cod/vp_pen_futuras.R")

Cálculo del VP de beneficios de pensiones futuras

t <- proc.time()
datos <- vp_pen_futuras(int_real, 0.02) 
proc.time() - t

##    user  system elapsed 
##   49.42    0.99   51.68

colSums(datos)/1e9

##          Inv          Pen          Suc          SEM          Cot       Cuotas 
##  22.30608415 144.85021090  33.26855290  15.75063412   0.04002766  59.58599631

Segundo modelo deterministico

Objetivo

Poder comparar el primer modelo (había salido mal)
calcular las pensiones por anualidades de mortalidad para más eficiencia

source("cod/anualidades.R")

Precálculo de los qx

pc_qx <- function(){
    w <- 115
  difsex <- list()
  for(sex in 1:2){
    tqx <- list()
    for(x in 20:w){
      res <- (w-x)
      qx <- sapply(1:res, function(y) all_qx[[sex]][x+y, y])
      tqx[[x-19]] <- qx
    }
    difsex[[sex]] <- tqx
  }
  return(difsex)
}

qx <- pc_qx()

Función para leer sumas de columnas de matrices triangulares

readt <- function(mat){
  return(sapply(1:ncol(mat), function(x) sum(mat[1:x,x])))
}

Cálculo con pensiones precalculadas

source("cod/vp_pen_fut.R")

t <- proc.time()
datos2 <- vp_pen_fut(int_real, 0.02)
proc.time() - t

##    user  system elapsed 
##    6.78    0.13    7.01

colSums(datos2)

##          Inv          Pen          Suc          SEM          Cot       Cuotas 
##  23882659791 137587976257  27864304544  14879120060     37473445  67691775052

Modelo estocástico

Precálculos

Obteniendo la probabilidad de meses para binomial

curva_densidad$media <- curva_densidad$media/12
curva_densidad$std <- curva_densidad$std/12

Valor presente de cotizaciones

cotizaciones <- cotizaciones*niveles

Cuotas de los cotizantes

cuotas <- list()
for(i in 1:5196){
  cuotas[[i]] <- sapply(27:30, function(x) sum(cotizaciones[i,1:12+(12*(x-1))] > 1e4))
}

Simplificando también la probabilidad de los hijos

oqx <- sapply(0:24, function(x)
  sapply(x:24, function(y)
    (all_qx[[1]][y+1,y-x+1] + all_qx[[2]][y+1,y-x+1])/2
    )
  )

Función estocástica

source("cod/proy_est.R")

system.time(est <- proy_est(int_real, 0.02))

##    user  system elapsed 
##   10.13    0.06   10.39

colSums(est)/1e9

##        Inv        Pen        Suc        SEM       InvF       PenF       SucF 
##   4.361823  17.350398   4.439264   2.055147  21.910805 146.840796  33.781063 
##       SEMF     CotAct     CotPen 
##  15.916279  68.543529   0.000000

Paralelización

variables <- c("cotizantes", "cotizaciones", "qx", "all_qix", "int_real", "proy_est", "cot_min_inv", "cuotas", "curva_densidad", "curva_salarial", "oqx", "pensionados")
proy_est_par <- function(n, cores) {
  # Crear un clúster seguro
  cl <- makeCluster(cores) #min(detectCores()/2,
  # Exportar las variables necesarias al clúster
  clusterExport(
    cl,
    varlist = variables,
    envir = environment()
  )
  # Ejecutar la proyección en paralelo
  resultados <- parSapply(cl, 1:n, function(x) {
    proy_est(int_real, 0.02)
  })
  # Detener el clúster después de la ejecución
  stopCluster(cl)
  return(resultados)
}

Guardado de las proyecciones

n <- 1000
# system.time(proyecciones <- proy_est_par(n, 4))
# copia <- lapply(1:n, function(x) as.data.frame(proyecciones[,x]))
# write.csv(copia, "res/datos.csv")

proyecciones <- read.csv("res/datos.csv")
proyecciones <- proyecciones %>% select(-1)
proyecciones <- lapply(1:1000, function(x) proyecciones[,1:10 + 10*(x-1)])

Media estocástica

sumas <- sapply(1:n, function(x) colSums(as.data.frame(proyecciones[[x]])))
(media <- sapply(1:10, function(y) mean(sumas[y,])))

##  [1] 4.314899e+09 1.702378e+10 4.513780e+09 2.031647e+09 2.279332e+10
##  [6] 1.463577e+11 3.326937e+10 1.590746e+10 7.047005e+10 4.647387e+04

Intervalo de confianza 99%

sapply(1:10, function(y) quantile(sumas[y,],0.99))/1e9

##        99%        99%        99%        99%        99%        99%        99% 
##   4.605287  18.009331   5.092541   2.077267  26.003750 149.350477  35.147569 
##        99%        99%        99% 
##  16.086308  98.375602   0.000000

sapply(1:10, function(y) quantile(sumas[y,],0.01))/1e9

##         1%         1%         1%         1%         1%         1%         1% 
##   3.974783  16.046940   3.959611   1.982637  19.949642 143.425566  31.202187 
##         1%         1%         1% 
##  15.730189  58.810836   0.000000

Convergencia a 3 cifras

sum((media*1000 + colSums(est)/1e9)/1001/media)/10

## [1] 0.999001

Gráfico de convergencia

Cálculo de pensiones actuales, determinístico

# IMPORTANTE: ninguna pensión sobrepasa los 2 millones
pensiones_actuales <- function(int, inf){
  v <- (1+inf)/(1+int)
  anual <- (1-v)/(v^(-1/12)-1)
  anmag <- (anual + v^(11/12)) 
  res <- data.frame(Inv = rep(0, 96), Pen = 0, Suc = 0, SEM = 0)
  for(id in 1:365){
    x <- pensionados$edad[id]
    sexo <- pensionados$sexo[id]
    tipo <- pensionados$COD_TIPO_PENSION[id]
    vp_pen <- pensionados$MONTO[id]*anmag
    if(tipo == "Invalidez" | tipo == "Vejez"){ # pension regular

      aqx <- qx[[sexo]][[x-19]]
      tpx <- c(1,cumprod(1-aqx))
      aqx <- c(aqx,1)
      if(tipo == "Invalidez" ){
        res$Inv <- res$Inv + c(tpx*vp_pen,rep(0, 95-(115-x)))
      } else {
        res$Pen <- res$Pen + c(tpx*vp_pen,rep(0, 95-(115-x)))
      }
      # viudez y orfandad

      res$Suc <- res$Suc + c(readt(viu[[3-sexo]][[x-19]]*tpx*aqx)*vp_pen,rep(0, 95-(115-x)))
      if(x < 50){

        res$Suc <- res$Suc + c(readt(orf[[x-19]]*tpx*aqx)*vp_pen,rep(0, 95-(115-x)))        
      }

    } else {
      if(pensionados$COD_PARENTESCO[id] == "H"){
        if(x >= 25){ # ya no más pensión
          next
        }
        res$Suc <- res$Suc + c(cumprod(1-oqx[[x+1]])*vp_pen, rep(0, 96-(25-x)))
      } else {
        res$Suc <- res$Suc + c(cumprod(1-qx[[sexo]][[x-20]])*vp_pen, rep(0, 95-(115-x)))
      }
    }
  }
  res$SEM <- (res$Inv + res$Pen + res$Suc)*anual/(anual+v^(11/12))*0.085
  return(res*v^(1:95))
}

pen_act <- pensiones_actuales(int_real, 0.02)

colSums(pen_act)/1e9

##       Inv       Pen       Suc       SEM 
##  4.270028 16.204969  4.840757  1.784558

Cálculo beneficio devengado, determinístico

Es el mismo código solo que añadiendo n_cot/cuota_ini

source("cod/ben_dev.R")

system.time(ben <- ben_dev(int_real, 0.02))

##    user  system elapsed 
##   54.64    1.01   56.78

colSums(ben)

##         Inv         Pen         Suc         SEM         Cot      Cuotas 
## 12633021826 77285552770 18560289404  8524945456    18099101 59585996314

Proyecciones demográficas, determinístico

Cotizantes actuales

source("cod/proy_dem.R")

demo <- proy_dem()
demo <- cbind(Año = 2025:2120, demo)

Recordando que después de 2054 no pueden haber hijos ya que se cumple la edad máxima Si hoy una persona tiene 21 años (edad mínima del fondo) entonces al 2054 tendrá 50 años, por lo que su hijo tiene 25 años y ya no recibe pensión

Pensionados actuales

demo_pen <- pen_dem()
demo_pen <- cbind(Año = 2025:2120, demo_pen)

Proyecciones financieras

datos <- as.data.frame(datos[c(2:96,1),])
rownames(datos) <- 1:96

Egresos

egresos <- data.frame(
  Año = 2025:2120,
  Curso = rowSums(pen_act[,1:3]),
  Futuras = rowSums(datos[,1:3]),
  SEM = rowSums(data.frame(pen_act$SEM, datos2$SEM))
) %>% mutate(
  Total = Curso + Futuras + SEM
)
write_xlsx(egresos, "res/egresos.xlsx")

Ingresos

ingresos <- data.frame(
  Año = 2025:2120,
  Cuotas = datos2$Cuotas,
  Cotizaciones = datos2$Cot,
  Reserva = 0,
  Inversiones = 0,
  Total = 0
)
ingresos$Reserva[1] <- 437031071554.9
ingresos$Inversiones[1] <- (ingresos$Reserva[1] + 
  (ingresos$Cuotas[1] + ingresos$Cotizaciones[1]-egresos$Total[1])/2)*int_real
ingresos$Total[1] <- ingresos$Cuotas[1] + ingresos$Cotizaciones[1] + ingresos$Inversiones[1]
for(i in 2:96){
  ingresos$Reserva[i] <- ingresos$Reserva[i-1] + ingresos$Total[i-1] - egresos$Total[i-1]
  ingresos$Inversiones[i] <- (ingresos$Reserva[i] + 
  (ingresos$Cuotas[i] + ingresos$Cotizaciones[i]-egresos$Total[i])/2)*int_real
  ingresos$Total[i] <- ingresos$Cuotas[i] + ingresos$Cotizaciones[i] + ingresos$Inversiones[i]
}

Reservas

reserva <- data.frame(Año = 2025:2120,
                      Salarios = datos2$Cuotas/0.15,
                      Reserva_Inicial = ingresos$Reserva,
                      Ingresos = ingresos$Total,
                      Egresos = egresos$Total,
                      Reserva_final = c(ingresos$Reserva[-1],
  ingresos$Reserva[96] + ingresos$Total[96] - egresos$Total[96])) %>% 
  mutate(
    Tasa_Ingresos = Ingresos/Salarios,
    Tasa_Costos = Egresos/Salarios
  )
ingresos <- ingresos %>% select(-Reserva)
write_xlsx(ingresos, "res/ingresos.xlsx")
write_xlsx(reserva, "res/reserva.xlsx")

Gráficos

Costo de pensiones

costo <- data.frame(
  Año = 2025:2120,
  Inv = rowSums(data.frame(datos2$Inv, pen_act$Inv))/1e6,
  Suc = rowSums(data.frame(datos2$Suc, pen_act$Suc))/1e6,
  Pen = rowSums(data.frame(datos2$Pen, pen_act$Pen))/1e6
)

Ingresos

Análisis de sensibilidad

Pesimista

colSums(vp_pen_futuras(int_real-0.005, 0.025))

##          Inv          Pen          Suc          SEM          Cot       Cuotas 
##  28065001059 189796240721  45799717612  20714458177     52654763  65501124289

colSums(pensiones_actuales(int_real-0.005, 0.025))

##         Inv         Pen         Suc         SEM 
##  4835362573 18082847978  7217618648  2169494435

Optimista

colSums(vp_pen_futuras(int_real+0.005, 0.015))

##          Inv          Pen          Suc          SEM          Cot       Cuotas 
##  17965524904 111921968445  24642632754  12147257854     30734474  54439380313

colSums(pensiones_actuales(int_real+0.005, 0.015))

##         Inv         Pen         Suc         SEM 
##  3802861300 14594739793  5067054024  1627283898

Matriz de sensibilidad

Razones de solvencia para 25 casos

act <- 437031071554.9+40930473298.54
# rs <- matrix(0,5,5)
# for(i in 1:5){
#   for(j in 1:5){
#     cot <- colSums(vp_pen_futuras(int_real+0.02*(i-3), 0.02+0.01*(j-3)))
#     pens <- colSums(pensiones_actuales(int_real+0.02*(i-3), 0.02+0.01*(j-3)))
#     rs[i,j] <- (act + sum(cot[5:6]))/(sum(cot[1:4]) + sum(pens))
#   }
# }

# razones <- as.data.frame(t(rs))
# razones <- rbind((int_real+0.02*(1:5-3))*100, razones)
# razones <- cbind(c(0,(0.02+0.01*(1:5-3))*100), razones)
# names(razones)[1] <- "Inf"
# write_xlsx(razones, "res/matriz_sens.xlsx")

read_excel("res/matriz_sens.xlsx", skip = 1)

## # A tibble: 5 × 6
##     `0` `3.793611978731801` `5.793611978731801` `7.793611978731802`
##   <dbl>               <dbl>               <dbl>               <dbl>
## 1     0               1.37                2.27                 3.53
## 2     1               1.03                1.76                 2.81
## 3     2               0.758               1.34                 2.21
## 4     3               0.554               1.01                 1.71
## 5     4               0.399               0.751                1.31
## # ℹ 2 more variables: `9.793611978731802` <dbl>, `11.7936119787318` <dbl>